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ABSTRACT 

Algorithms for the generation of three-dimensional unstructured surface and volume 
grids are discussed. These algorithms are based on incremental insertion and local opti- 
mization. The present algorithms are very general and permit local grid optimization based 
on various measures of grid quality. This is very important; unlike the two-dimensional 
Delaunay triangulation, the 3-D Delaunay triangulation appears not to have a lexico- 
graphic characterization of angularity. (The Delaunay triangulation is known to minimize 
that maximum containment sphere, but unfortunately this is not true lexicographically). 
Consequently, Delaunay triangulations in three-space can result in poorly shaped tetrahe- 
dral elements. Using the present algorithms, three-dimensional meshes can be constructed 
which optimize a certain angle measure, albeit locally. We also discuss the combinatorial 
aspects of the algorithm as well as implementational details. 
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INTRODUCTION 


Unstructured grids offer great advantages in the automatic generation of computa- 
tional grids about complex geometries. In the present work we consider triangulation 
methods in two and three space dimensions. In order to minimize the amount of re- 
quired user interaction, incremental triangulation algorithms are employed which permit 
automatic (adaptive) point placement. Several incremental algorithms exist for the trian- 
gulation of two and three dimensional domains. Most of these algorithms are designed to 
produce the Delaunay triangulation of specified sites which form the vertices of the mesh. 
In the present work, we have chosen an incremental triangulation algorithm based on point 
insertion and local edge swapping. The advantage of this algorithm is that the edge swap- 
ping can be used to optimize user specified grid quality measures. (This algorithm can 
be used to generate a Delaunay triangulation by proper choice of measure.) In addition, 
we produce Steiner triangulations by inserting new sites at triangle circumcenters of the 
existing triangulation to minimize the maximum triangle aspect ratio (or any other user 
specified criteria). An objective of the current research is to extend this two-dimensional 
algorithm to include the generation of surface grids or volume grids in three-dimensions. 
The algorithm based on edge swapping is well known in two dimensions [1] but the exten- 
sion to surface grids and volume grids in three dimensions in not straightforward. In the 
present work, we show the extension of edge swapping to include the generation of surface 
grids and three dimensional volume grids. Differing edge swapping optimization criterion 
are considered. Steiner triangulations in three dimensions are produced by placing new 
sites at circumsphere centers of the existing triangulation. 


Motivation: 

• Develop fast automatic three-dimensional unstructured grid generation capability for 
complex geometries. 

Objectives: 

• Generalize a 2-D unstructured grid generation algorithm based on incremental inser- 
tion and edge swapping. 

• Generation of surface grids on 2-manifolds in 3-D. 

• Automatic adaptation to maintain surface fidelity in 3-D. 

• Extend incremental insertion and edge swapping to 3-D volume mesh generation. 

• Investigation of differing edge swapping optimization criterion. 

• Automatic point insertion (Steiner triangulations). 
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Edge Swapping or Watson’s Algorithm? 


Conceptually speaking, the incremental algorithm based on edge swapping differs 
significantly from the incremental Watson algorithm [4] for the Delaunay triangulation. 
However, when edge swapping is used to generate a Delaunay triangulation, both methods 
have similar implementations. Both methods begin by inserting a new site P in the interior 
of an existing triangulation. The edge swapping algorithm begins by connecting site P 
to the existing triangulation. The union of all triangles incident to P forms a polygon 
surrounding this new site. The edges of this polygon must each pass the circumsphere 
test. Edges failing this test must be swapped. The process continues until all edges of 
the polygon pass. The algorithm is trivial to implement using recursive programming. 
Conceptually, Watson’s algorithm finds all triangles with circumcircles containing site P. 
This identifies invalid edges that when removed reveal a polygon with vertices visible from 
P. The vertices of this polygon are then connected to P. The actual implementation of 
Watson’s algorithm is similar to the edge swapping code. The first step is to find any 
triangle with its circumcircle containing site P. Once this triangle is found, edges are 
reconnected and neighboring triangles tested. This can also be easily programmed 
recursively. The resulting code is very similar to the edge swapping implementation. 

Edge Swapping Watson’s Algorithm 



Recursive Edge Swapping on Suspect Edges Remove Invalid Edges and Reconnect 
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Min-Max and Max-Min Triangulations 


The edge swapping algorithm easily extends to allow edge swapping based on criteria 
other than the circumsphere test. The circumsphere test is known to be equivalent to the 
angularity optimization that chooses the diagonal position that maximizes the minimum 
interior angle (Max-Min triangulation). Another popular method chooses a diagonal po- 
sition for a given convex triangle pair that minimizes the maximum interior angle (the 
Min-Max triangulation). If the diagonal position is swapped using this measure then the 
four neighboring edges must be checked (the Max-Min version of this algorithm only re- 
quires that two edges be checked). It should also be noted that the Min-Max triangulation 
can produce triangulations that locally minimize the maximum angle but are not globally 
optimal with respect to the same measure. 


• Advantage of incremental edge swapping is that it permits other types of local opti- 
mization. 

• Requires more sophisticated recursion implementation and usually only produces 
a locally optimum triangulation. 

• Typical local optimizations minimize the maximum interior angle (Min-Max) or 
maximize the minimum interior angle (Max-Min, Delaunay triangulation). 



Delaunay (Max-Min) Configuration 


Min-Max Configuration 


• Min-Max triangulation favors less obtuse triangles (important for highly stretched 
meshes). 
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Steiner Triangulations 


A Steiner triangulation is any triangulation that adds additional sites into an existing 
triangulation to improve a specified measure of grid quality. For the present application, 
Delaunay triangulations are constructed which refine a given triangle by adding a new 
site at the circumcenter of the triangle. This technique is well known (Homes and Snyder 
[3], Warren et al [4]). The choice of circumcenter guarantees that no other point in the 
triangulation lies closer than the radius of this circle. The initial triangulation consists of 
a Delaunay triangulation of the boundary points (possibly edge swapped to maintain the 
boundary curve definitions). Points are inserted to minimize the maximum triangle aspect 
ratio using a dynamic heap data structure. 


• Steiner Triangulation - insert additional sites into an existing triangulation to improve 
a specified measure of grid quality. 



Initial Triangulation After Steiner Point Insertion 

and Interior Edge Removal 

• Insert sites at circumcenters to minimize triangle aspect ratio. 
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3-D Surface Patch Projection 


Three-dimensional surface grids are constructed from rectangular surface patches (as- 
sumed at least C° smooth) using a generalization of the 2-D Steiner triangulation scheme. 
Points are placed on the perimeter of each patch using an adaptive refinement strategy 
based on absolute error and curvature measures. The surface patches are projected onto 
the plane. Simple stretching of the rectangular patches permits the user to produce pref- 
erentially stretched meshes. (This is useful near the leading edge of a wing for example.) 


• Sites are distributed on the perimeter of each patch using 1-D adaptive subdivision. 

• Curve adaptation criterion: 

• Absolute tolerance 

• Curvature tolerance 

• Patches projected into 2-space for triangulation. 

• Mapping to rectangle permits simple stretching which results in directional bias- 
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Adaptive Patch Triangulation 


The triangulation takes place in the two dimensional ( s , t ) plane. The triangulation 
is adaptively refined using Steiner point insertion to minimize the maximum user specified 
absolute error and curvature tolerance on each patch. The absolute error is approximated 
by the perpendicular distance from the triangle centroid (projected back to 3-space) to the 
true surface. The user can further refine based on triangle aspect ratio in the ( s,t ) plane 
if desired. 


• Adaptive Delaunay triangulation in ( s,t ) plane. 

• Adaptive insertion based on absolute error and maximum curvature tolerance. 



• Calculation of triangulation absolute error by measurement of distance from face 
centroid to true surface. 

• Optional refinement based on triangle aspect ratio. 
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Boeing 737 with Flaps Deployed 
Adaptive Surface Triangulation 



3-D Incremental Insertion and Optimization 

The 3-D edge swapping methodology is based on Lawson’s result [5] that there exists 
at most two ways to triangulate n + 1 points in R n . It can be shown that if a configuration 
of 5 points in 3-D allows only one triangulation, it will satisfy the Delaunay circumsphere 
test. If a configuration of 5 points in 3-D allows two triangulations, then only one will 
satisfy the Delaunay circumsphere test. The edge swapping algorithm is based on locally 
swapping convex non-Delaunay triangulations of 5 points into its Delaunay counterpart. 
It has been shown that during the optimization process for a general mesh the procedure 
may get stuck in a local optimum. However, Rajan [1] proved that inserting a site into 
an existing Delaunay triangulation followed by edge swapping is guaranteed to recover the 
new Delaunay triangulation. This forms the basis of our 3-D Incremental Insertion and 
Optimization algorithm for 3-D Volume Mesh Generation. 

• Combinatorial theorem (Lawson): There exist at most two ways to triangulate n + 1 
points in R n (subject to convexity). 



2 Tetrahedron Configuration 3 Tetrahedron Configuration 


• Edge swapping an arbitrary 3-D triangulation (subject to circumsphere optimization) 
is not sufficient to guarantee that a Delaunay triangulation will be produced. 

• Site insertion into an existing Delaunay triangulation followed by edge swapping will 
produce a new Delaunay triangulation (Rajan, 1991). 
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Local and Global Optimization 


As long as a metric can be defined for the two triangulations of a configuration of 5 
points in three dimensions, edge swapping locally to optimize that metric can be performed. 
In this way, edge swapping provides a general framework which allows for different criteria 
of judging mesh quality. This permits criteria other than the Delaunay circumsphere 
test. For instance, we can maximize the minimum face (or solid) angle, or minimize the 
maximum face (or solid) angle, or reduce the number of edges and volumes in the mesh. 
Local optimization is not guaranteed to yield the globally optimal mesh, but typically it 
produces significantly improved meshes. 


• 3-D edge swapping permits local optimization based on many possible criterion such 
as: 

• Circumsphere 

• Maximize the minimum face (or solid) angle 

• Minimize the maximum face (or solid) angle 

• Minimize the total number of edges and volume 

• Local optimization usually does not produce a globally optimal triangulation but can 
significantly improve the mesh quality. 



Histogram showing improvement in maximum face angle using edge swapping procedure 
on triangulation of 500 random sites in unit cube. 
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3-D Volume Triangulation 


In order to perform our 3-D Incremental Insertion and Optimization algorithm for 3-D 
Volume Mesh Generation, we need a constrained Delaunay triangulation of the boundary 
in which we can incrementally add new sites. This Constrained Delaunay triangulation 
of the boundary sites is constructed using the Tanamura algorithm. Steiner points are 
introduced into this triangulation in a way that optimizes some criterion of mesh quality. 
We have been working with trying to minimize the maximum cell (tetrahedron) aspect 
ratio by introducing a site at the circumcenter of the tetrahedron. It should be noted 
that construction of constrained Delaunay triangulation is not always possible and we may 
need to add Steiner points to the interior of the untetrahedralizable volume to make it 
tetrahedralizable. 


• Constrained Delaunay triangulation of boundary sites is constructed using Tanamura 
algorithm implemented by Marshal Merriam at NASA Ames 

• Steiner points are inserted into this triangulation which minimize the maximum tetra- 
hedron aspect ratio (or any other user specified criterion). 

• Caveat Emptor: Construction of constrained Delaunay triangulation of boundary sites 
is not always possible unless Steiner points are added to the interior of the volume. 
Example: 



Nontetrahedralizable twisted prism with reflex exterior faces 
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CONCLUSIONS 


We have demonstrated the feasibility of adaptive Steiner point insertion and edge 
swapping to generate surface meshes about complex geometries. Edge swapping algorithms 
have the advantage over certain other standard triangulation techniques because of the 
ability to improve mesh quality through local optimization. We have also shown how 
to extend the edge swapping idea to three dimensions using the combinatorial result of 
Lawson. The resulting algorithm permits adaptive Steiner triangulations about complex 
geometries. 


• Demonstrated the feasibility of adaptive Steiner point insertion and edge swapping to 
generate surface meshes about complex geometries. 

• Implemented an incremental insertion and edge swapping algorithm to generate 3-D 
volume meshes. 

• Demonstrated the unique capability of edge swapping to improve mesh quality using 
various optimization criteria. 


REFERENCES 

1. Green, P.J. and Sibson, R., “Computing the Dirichlet Tesselation in the Plane”, The 
Computer Journal, Vol. 21, No. 2, 1977, pp. 168-173. 

2. Rajan, V.T. Optimality of the Delaunay Triangulation in R d . In Proceedings of the 
7th Annual ACM Symposium on Computational Geometry , 1991, pp. 357-372. 

3. Holmes, G. and Snyder, D., “The Generation of Unstructured Triangular Meshes 
using Delaunay Triangulation,” in Numerical Grid Generation in CFD , pp. 643-652, 
Pineridge Press, 1988. 

4. Warren, G., Anderson, W.K., Thomas, J.L., and Krist, S.L.,“Grid Convergence for 
Adaptive Methods,”, AIAA paper 91-1592-CP, Honolulu, Hawaii, June 24-27,1991. 

5. Lawson, C. L., “Properties of n-dimensional Triangulations” CAGD, Vol. 3, April 
1986, pp. 231-246. 


461 



